function calculateSimilarityMaps(Params,nrBand,nrSession)

% This function calculates time-varying and across-session
% intersubject similarity maps. Results can be accessed
% through memMaps.resultMap -field.
%
% Inputs:
% Params - struct containing all necessary parameters
% nrBand - frequency subband index (note: 0 refers to full band)
% nrSession - session index
%
% Example:
% load analysisParameters
% calculateSimilarityMaps(analysisParameters,4,1);
% Then access ISC-map of 10th time-interval:
% load memMaps
% ISCmap10 = memMaps.resultMap.win.band4.Session1.cor.Data(10).xyz;
%
% See also:
% INITPARAMS
% RUNANALYSIS


% Last modified: 10.11.2010 by Jukka-Pekka Kauppi
% Tampere University of Technology
% Department of Signal Processing
% e-mail: jukka-pekka.kauppi@tut.fi

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


showTime(1);

Pub = Params.PublicParams;
Priv = Params.PrivateParams;


% check parameters:

if ~Pub.calcStandard
    disp('Standard ISC analysis not specified, canceling computation of similarity maps...')
    return
end

load([Pub.dataDestination 'memMaps'])

if Pub.nrFreqBands > 0
    mMapFilt = memMaps.(Priv.filtMapName).([Priv.prefixSession num2str(nrSession)]);
end

mMapOrig = memMaps.(Priv.origMapName).([Priv.prefixSession num2str(nrSession)]);
mMapResult = memMaps.(Priv.resultMapName);

try
if ~isempty(mMapResult)
    if isfield(mMapResult,'whole')
        if isempty(mMapResult.whole) == false
            if Pub.corOn
                if mMapResult.whole.([Priv.prefixFreqBand num2str(nrBand)]).(...
                        [Priv.prefixSession num2str(nrSession)]).cor.Writable == false
                    disp('Similarity maps written already, canceling computation...')
                    return
                end
            end
            if Pub.ssiOn
                if mMapResult.whole.([Priv.prefixFreqBand num2str(nrBand)]).(...
                        [Priv.prefixSession num2str(nrSession)]).ssi.Writable == false
                    disp('Similarity maps written already, canceling computation...')
                    return
                end
            end
            if Pub.nmiOn
                if mMapResult.whole.([Priv.prefixFreqBand num2str(nrBand)]).(...
                        [Priv.prefixSession num2str(nrSession)]).nmi.Writable == false
                    disp('Similarity maps written already, canceling computation...')
                    return
                end
            end
            if Pub.kenOn
                if mMapResult.whole.([Priv.prefixFreqBand num2str(nrBand)]).(...
                        [Priv.prefixSession num2str(nrSession)]).ken.Writable == false
                    disp('Similarity maps written already, canceling computation...')
                    return
                end
            end
        end
    end
    if isfield(mMapResult,'win')
        if isempty(mMapResult.win) == false
            if Pub.corOn
                if mMapResult.win.([Priv.prefixFreqBand num2str(nrBand)]).(...
                        [Priv.prefixSession num2str(nrSession)]).cor.Writable == false
                    disp('Similarity maps written already, canceling computation...')
                    return
                end
            end
            if Pub.ssiOn
                if mMapResult.win.([Priv.prefixFreqBand num2str(nrBand)]).(...
                        [Priv.prefixSession num2str(nrSession)]).ssi.Writable == false
                    disp('Similarity maps written already, canceling computation...')
                    return
                end
            end
            if Pub.nmiOn
                if mMapResult.win.([Priv.prefixFreqBand num2str(nrBand)]).(...
                        [Priv.prefixSession num2str(nrSession)]).nmi.Writable == false
                    disp('Similarity maps written already, canceling computation...')
                    return
                end
            end
            if Pub.kenOn
                if mMapResult.win.([Priv.prefixFreqBand num2str(nrBand)]).(...
                        [Priv.prefixSession num2str(nrSession)]).ken.Writable == false
                    disp('Similarity maps written already, canceling computation...')
                    return
                end
            end
        end
    end
else
     disp('Memory pointers not available, calculation canceled...')
    return
end
catch
   error('Cannot update similarity maps possibly due to changes in analysis parameters. To update the data, delete bin-files in results-folder and try analysis again!')
   return
end


clear memMaps

[pl,ms,en] = computer;

cDat = zeros([Priv.dataSize(nrSession,[4 2 3]), Priv.nrSubjects]);

% load mask:
if strcmp(Pub.fileFormatSubj,'nii')
    bmask = load_nii(Priv.brainMask);
    bmask = single(bmask.img);
elseif strcmp(Pub.fileFormatSubj,'mat')
    bmask = load(Priv.brainMask);
    fiel = fields(bmask);
    bmask = bmask.(fiel{1});
    bmask = single(bmask);
else
    error('Mask must be mat- or nii-file!')
    return
end
bmask = logical(bmask);

INDS = [];
for hh = 1:Priv.nrSubjects-1
    INDS = [INDS (1+hh*Priv.nrSubjects):(...
        1+hh*Priv.nrSubjects+hh-1)];
end
iter = 0;

if Priv.nrTimeIntervals(nrSession) == 0
    nrFrames = 1;
else
    nrFrames = Priv.nrTimeIntervals(nrSession);
end



%%%%%%%%%%%%%%%%%%%%
%warning('Modified window steps!!!!')
%clear sI;for k=1:24;if mod(k,5)==0;sI(k)=8;else;sI(k)=9;end;end;sI=cumsum(sI)-8
%Priv.startInds{1} = sI;
%clear sI;for k=1:39;if mod(k,5)==0;sI(k)=8;else;sI(k)=9;end;end;sI=cumsum(sI)-8
%Priv.startInds{1} = sI;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calculation of similarity maps:

tic
for xx = 1:Priv.dataSize(nrSession,1)
    disp(['x:' num2str(xx) '/' num2str(Priv.dataSize(nrSession,1))])
    % process only non-zero slices:
    if sum(sum(squeeze(bmask(xx,:,:)))) > 0
        cDat = zeros(Priv.dataSize(nrSession,4),Priv.dataSize(nrSession,2),...
            Priv.dataSize(nrSession,3),Priv.nrSubjects);
        
        % get mapped source data of the subjects:
        for k = 1:Priv.nrSubjects
            if nrBand == 0 % load full-band data
                if ~strcmp(Priv.computerInfo.endian,en)
                    cDat(:,:,:,k) = swapbytes(mMapOrig.([Priv.prefixSubject num2str(k)]).Data(xx).tyz);
                else
                    cDat(:,:,:,k) = mMapOrig.([Priv.prefixSubject num2str(k)]).Data(xx).tyz;
                end
            else % load sub-band data
                if ~strcmp(Priv.computerInfo.endian,en)
                    cDat(:,:,:,k) = swapbytes(mMapFilt.([Priv.prefixSubjectFilt ...
                        num2str(k)]).([Priv.prefixFreqBand num2str(nrBand)]).Data(xx).tyz);
                else
                    cDat(:,:,:,k) = mMapFilt.([Priv.prefixSubjectFilt ...
                        num2str(k)]).([Priv.prefixFreqBand num2str(nrBand)]).Data(xx).tyz(:,:,:);
                end
            end
        end
        
        % init temp data matrices:
        kenDatawin = zeros([Priv.dataSize(nrSession,2:3),Priv.nrTimeIntervals(nrSession)]);
        corDatawin = zeros([Priv.dataSize(nrSession,2:3),Priv.nrTimeIntervals(nrSession)]);
        ssiDatawin = zeros([Priv.dataSize(nrSession,2:3),Priv.nrTimeIntervals(nrSession)]);
        nmiDatawin = zeros([Priv.dataSize(nrSession,2:3),Priv.nrTimeIntervals(nrSession)]);
        kenData = zeros([Priv.dataSize(nrSession,2:3)]);
        corData = zeros([Priv.dataSize(nrSession,2:3)]);
        ssiData = zeros([Priv.dataSize(nrSession,2:3)]);
        nmiData = zeros([Priv.dataSize(nrSession,2:3)]);
        for yy = 1:Priv.dataSize(nrSession,2)
            for zz = 1:Priv.dataSize(nrSession,3)
                if bmask(xx,yy,zz)
                    iter = iter + 1;
                    if mod(iter,10000) == 0
                        disp(['iter: ' num2str(iter) '/' num2str(length(find(bmask)))])
                        toc
                        tic
                    end
                    % obtain each subject's time series:
                    ts = squeeze(cDat(:,yy,zz,:));
                    % ikkunointi:
                    for wfr = 1:nrFrames
                        % calculate across-whole-session similarity values:
                        if wfr == 1
                            N = Priv.dataSize(nrSession,4);
                            %tic
                            if Pub.ssiOn
                                % Sign-similarity index:
                                dts = sign(diff(ts));
                                ssiData(yy,zz) = single(sum(abs(sum(dts')))/prod(size(dts)));
                            end
                            if Pub.nmiOn
                                % Mutual information of binary data:
                                dts(dts == -1) = 0;
                                % call binary_mi.m:
                                nmiData(yy,zz) = single(binary_mi(dts));
                            end
                            if Pub.corOn
                                % correlation coefficient calculation (ref: Matlab corrcoef.m)
                                [n1,m1] = size(ts);
                                xc = ts - repmat(sum(ts)/n1,n1,1);  % Remove mean
                                c1 = (xc' * xc) / (n1-1); % calculate inner products
                                d1 = sqrt(diag(c1)); % sqrt first to avoid under/overflow
                                dd = d1*d1';
                                dd(1:m1+1:end) = diag(c1); % remove roundoff on diag
                                r1 = c1 ./ dd;
                                r1 = r1(INDS);
                                % save all correlations:
                                corMatData(yy,zz,:) = r1;
                                % save mean correlation:
                                corData(yy,zz) = single(mean(r1));
                                
                                if mod(iter,10000) == 0
                                    disp(['Intersubject correlation value ('...
                                        num2str(xx) ',' num2str(yy) ',' ...
                                        num2str(zz) '): ' num2str(corData(yy,zz))])
                                end
                            end
                            % Kendall's coefficient of concordance:
                            % tic
                            if Pub.kenOn
                                [YY,R] = sort(ts);
                                RS = sum(R,2);
                                S = sum(RS.^2)-N*mean(RS).^2;
                                F = Priv.nrSubjects*...
                                    Priv.nrSubjects*(N*N*N-N);
                                kenData(yy,zz) = single(12*S/F);
                                %kenT = toc
                            end
                        end
                        
                        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                        % calculate similarity values for time-frames:
                        
                        if Priv.nrTimeIntervals(nrSession) > 0
                            N = Pub.windowSize;
                            ts_win = ts(Priv.startInds{nrSession}...
                                (wfr):Priv.startInds{nrSession}...
                                (wfr)+Pub.windowSize-1,:);
                            if Pub.ssiOn
                                % Sign-similarity index:
                                dts = sign(diff(ts_win));
                                ssiDatawin(yy,zz,wfr) = single(sum(abs(sum(dts')))/prod(size(dts)));
                            end
                            if Pub.nmiOn
                                % Mutual information of binary data:
                                dts(dts == -1) = 0;
                                nmiDatawin(yy,zz,wfr) = single(binary_mi(dts));
                            end
                            if Pub.corOn
                                % Mean of pairwise correlation (ref: Matlab corrcoef.m):
                                [n1,m1] = size(ts_win);
                                xc = ts_win - repmat(sum(ts_win)/n1,n1,1);  % Remove mean
                                c1 = (xc' * xc) / (n1-1); % calculate inner products
                                d1 = sqrt(diag(c1)); % sqrt first to avoid under/overflow
                                dd = d1*d1';
                                dd(1:m1+1:end) = diag(c1); % remove roundoff on diag
                                r1 = c1 ./ dd;
                                r1 = r1(INDS);
                                % save all correlations:
                                corMatDatawin(yy,zz,:,wfr) = r1;
                                % save mean correlation:
                                corDatawin(yy,zz,wfr) = single(mean(r1));
                            end
                            if Pub.kenOn
                                % Kendall's coefficient of concordance:
                                [YY,R] = sort(ts_win);
                                RS = sum(R,2);
                                S = sum(RS.^2)-N*mean(RS).^2;
                                F = Priv.nrSubjects*...
                                    Priv.nrSubjects*(N*N*N-N);
                                kenDatawin(yy,zz,wfr) = single(12*S/F);
                            end
                        end
                    end
                end
            end
        end
        
        % save results:
        %    tic % result saving time may vary considerably
        if Pub.ssiOn
            mMapResult.whole.([Priv.prefixFreqBand...
                num2str(nrBand)]).([Priv.prefixSession num2str(nrSession)]...
                ).ssi.Data.xyz(xx,:,:) = ssiData;
        end
        if Pub.nmiOn
            mMapResult.whole.([Priv.prefixFreqBand...
                num2str(nrBand)]).([Priv.prefixSession...
                num2str(nrSession)]).nmi.Data.xyz(xx,:,:) = nmiData;
        end
        if Pub.corOn
            mMapResult.whole.([Priv.prefixFreqBand...
                num2str(nrBand)]).([Priv.prefixSession...
                num2str(nrSession)]).cor.Data.xyz(xx,:,:) = corData;
        end
        if Pub.kenOn
            mMapResult.whole.([Priv.prefixFreqBand...
                num2str(nrBand)]).([Priv.prefixSession ...
                num2str(nrSession)]).ken.Data.xyz(xx,:,:)= kenData;
        end
        if nrFrames > 1
            for wfr = 1:size(ssiDatawin,3)
                if Pub.ssiOn
                    mMapResult.win.([Priv.prefixFreqBand...
                        num2str(nrBand)]).([Priv.prefixSession...
                        num2str(nrSession)]).ssi.Data(wfr).xyz(xx,:,:) = ssiDatawin(:,:,wfr);
                end
                if Pub.nmiOn
                    mMapResult.win.([Priv.prefixFreqBand...
                        num2str(nrBand)]).([Priv.prefixSession...
                        num2str(nrSession)]).nmi.Data(wfr).xyz(xx,:,:) = nmiDatawin(:,:,wfr);
                end
                if Pub.corOn
                    mMapResult.win.([Priv.prefixFreqBand...
                        num2str(nrBand)]).([Priv.prefixSession...
                        num2str(nrSession)]).cor.Data(wfr).xyz(xx,:,:) = corDatawin(:,:,wfr);
                end
                if Pub.kenOn
                    mMapResult.win.([Priv.prefixFreqBand...
                        num2str(nrBand)]).([Priv.prefixSession...
                        num2str(nrSession)]).ken.Data(wfr).xyz(xx,:,:) = kenDatawin(:,:,wfr);
                end
            end
        end
        %      toc
    end
end

% set maps non-writable:
if Pub.corOn
    load([Pub.dataDestination 'memMaps'])
    memMaps.resultMap.whole.([Priv.prefixFreqBand...
        num2str(nrBand)]).([Priv.prefixSession...
        num2str(nrSession)]).cor.Writable = false;
    save([Pub.dataDestination 'memMaps'],'memMaps')
    clear memMaps
    if Priv.nrTimeIntervals(nrSession) > 0
        load([Pub.dataDestination 'memMaps'])
        memMaps.resultMap.win.([Priv.prefixFreqBand...
            num2str(nrBand)]).([Priv.prefixSession...
            num2str(nrSession)]).cor.Writable = false;
        save([Pub.dataDestination 'memMaps'],'memMaps')
        clear memMaps
    end
end
if Pub.kenOn
    load([Pub.dataDestination 'memMaps'])
    memMaps.resultMap.whole.([Priv.prefixFreqBand...
        num2str(nrBand)]).([Priv.prefixSession...
        num2str(nrSession)]).ken.Writable = false;
    save([Pub.dataDestination 'memMaps'],'memMaps')
    clear memMaps
    if Priv.nrTimeIntervals(nrSession) > 0
        load([Pub.dataDestination 'memMaps'])
        memMaps.resultMap.win.([Priv.prefixFreqBand...
            num2str(nrBand)]).([Priv.prefixSession...
            num2str(nrSession)]).ken.Writable = false;
        save([Pub.dataDestination 'memMaps'],'memMaps')
        clear memMaps
    end
end
if Pub.ssiOn
    load([Pub.dataDestination 'memMaps'])
    memMaps.resultMap.whole.([Priv.prefixFreqBand...
        num2str(nrBand)]).([Priv.prefixSession...
        num2str(nrSession)]).ssi.Writable = false;
    save([Pub.dataDestination 'memMaps'],'memMaps')
    clear memMaps
    if Priv.nrTimeIntervals(nrSession) > 0
        
        load([Pub.dataDestination 'memMaps'])
        memMaps.resultMap.win.([Priv.prefixFreqBand...
            num2str(nrBand)]).([Priv.prefixSession...
            num2str(nrSession)]).ssi.Writable = false;
        save([Pub.dataDestination 'memMaps'],'memMaps')
        clear memMaps
    end
end
if Pub.nmiOn
    load([Pub.dataDestination 'memMaps'])
    memMaps.resultMap.whole.([Priv.prefixFreqBand...
        num2str(nrBand)]).([Priv.prefixSession...
        num2str(nrSession)]).nmi.Writable = false;
    save([Pub.dataDestination 'memMaps'],'memMaps')
    clear memMaps
    
    if Priv.nrTimeIntervals(nrSession) > 0
        load([Pub.dataDestination 'memMaps'])
        memMaps.resultMap.win.([Priv.prefixFreqBand...
            num2str(nrBand)]).([Priv.prefixSession...
            num2str(nrSession)]).nmi.Writable = false;
        save([Pub.dataDestination 'memMaps'],'memMaps')
        clear memMaps
    end
end


showTime(0);
